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1. Introduction 

Differential equations with a time-delayed feedback have attracted a lot of attention in 
various fields in recent years [1-7]. Due to their non-Markovian dynamics such equations 
exhibit very interesting phenomena. However, our understanding of their mathematical 
properties is still at the beginning [8,9]. This applies not only to deterministic systems 
which are more frequently studied but also to stochastic differential equations with a 
time-delayed feedback. 

Most of the differential equations with a time-delayed feedback that have been 
discussed so far may be written in the form 

x{t) = F[x(t)] + G[x(t-r)], (1) 

where F and G are certain functions. Here the quantity x(t) may be understood as a one- 
dimensional stream of data that evolves according to a non-local dynamics. However, 
slicing the time line into equidistant intervals of width r and arranging these segments 
line by line in a two-dimensional plane, the time-delayed feedback becomes local as it 
connects points of two neighboring intervals at the same relative position. Combined 
with an appropriate continuum limit this means that a one-dimensional system with a 
time-delayed dynamics can be mapped onto a 1+1-dimensional system on a strip with 
local dynamics. This correspondence works also for stochastic systems, where F or G 
involve random variables. 
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Figure 1. Mapping a one-dimensional time series with a constant time-delayed 
feedback onto a 1+1-dimensional strip. The time line is sliced into equal 
intervals of width r which are arranged line by line on a strip with spiral-like 
boundary conditions. As can be seen, the non-local feedback (as examplifled 
by the red arrow) turns into a local interaction between neighboring segments. 



The purpose of this work is to show that this mapping can be used to relate 
certain one-dimensional time-delayed stochastic differential equations to well-known 
1+1-dimensional stochastic models with local dynamics and thereby predict their 
statistical properties. As an example we study the stochastic differential equation 

x(t) = ax{t) + bx(t - r) + cx(t) f (t) , (2) 

where x(t) ei is the one-dimensional time series, a, b and c are real- valued parameters, 
r > is the temporal delay and denotes a Gaussian white noise with the correlation 

mm) = s(t - v . (3) 

Recently this equation was studied numerically by Brettschneider [10], who observed 
extremely strong fluctuations of x(t). By relating this system to a 1+1-dimensional 
stochastic differential equation we will show that these fluctuations can be described in 
terms of the Kardar-Parisi-Zhang (KPZ) universality class [11]. 

Since is multiplied by the field x(t), the noise in Eq. (J2J) is said to be 
multiplicative. As in any system with multiplicative noise the iteration scheme has 
to be specified. Throughout this paper we shall assume that Eq. (J2J) is iterated in the 
Ito sense. 

Before starting we note that the stochastic differential equation @ has the following 
symmetry properties: 

(a) The equation is linear and invariant under translations in time. 

(b) A change of the time scale t — > At rescales the parameters by 

a — ► \a , b — > Xb , c— >^/\c, r — ► At (4) 

(c) Since the noise is statistically symmetric under reflections, the equation is invariant 
under a change of sign c — ► — c. 

(d) Under the transformation x(t) — > e^xit) the parameters change as 

a— >a — jj,, b^be~^ T (5) 
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In what follows we use the invariances (b) and (c) to set c = 1 without loss of generality. 
Moreover, the property (d) can be used to fix the parameter a. However, for the sake 
of clarity we keep a as a free parameter throughout the calculations. 

Since the stochastic differential equation invokes a temporally delayed feedback r, 
the initial state of the evolution is specified by the whole history of x(t) in a finite 
temporal interval of width r, for example t G [0, r). For simplicity we shall assume that 
x(t) = 1 on this interval. The question posed here is how x{t) is distributed for large t. 

Unless stated otherwise we assume that the parameter b is non-negative. This 
ensures that x(t) with a positive initial condition cannot change sign. To see this consider 
the transformation x(t) = e y ^ which leads to the stochastic differential equation 
y(t) = a + bert-^-v® + f (t) with a non-multiplicative noise and the initial condition 
y(t) = for < t < t. Obviously, for b > this differential equation cannot evolve 
towards y{t) = — oo within finite time, hence x(t) is positive for all t. 

The paper is organized as follows. In Sect. 2 the stochastic differential equation (J2J) 
is discretized and integrated numerically. In Sect. 3 we show how the one-dimensional 
iteration with non-local delay can be mapped onto a two-dimensional system with a 
local update rule. Carrying out the continuum limit this update rule turns out to be 
equivalent to a KPZ equation. Finally in Sect. 4 we compare the theoretical prediction 
with the numerical results. We finish the paper with some conclusions. 

2. Numerical analysis 

To study Eq. (j2J) numerically we discretize the time variable t = nh in steps of h. Then 
(using the Ito scheme) the discretized version of Eq. ((2j) reads 

x n = (l + ah)x n -. 1 + bhx n -. k + vhx n -iz n - 1 (6) 

where k — r/h is the discretized delay. The Z n S £1X6 independent Gaussian random 
numbers with variance 1 which may be generated by a Box-Miiller transformation [12]. 
Since a Gaussian distribution is unbound it can happen that the random number z n -\ 
takes on a very large negative value so that x n+ i becomes negative. This overshooting to 
negative values is a numerical artifact caused by the discretization since in the contiuum 
limit x(t) is always positive as long as b > 0. There are several ways to deal with this 
problem: 

(a) The time increment h may be decreased until the probability for overshooting to 
negative values is practically zero. However, this method is often inefficient and 
unsafe. 

(b) Another possiblity is to discard all those updates that would produce negative 
values. This method is safe but it introduces systematic numerical errors in form 
of an effective bias towards larger amplitudes of x n . 
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Figure 2. Distribution of y n = log(x n ) of the time series x n with delay k = 100 
after n = 10 5 updates for c = 1 and a = 1/2. For 6 = one obtains a 
Gaussian distribution centered at the origin. For 6 = 1, however, the center of 
the distribution moves to the right as time proceeds. Moreover, the distribution 
is no longer Gaussian, instead it becomes skew. 



(c) The most accurate method is a semi-analytical approach suggested by Dornic et al. 
in Ref. [13]. 

For the purpose of the present work, where we are primarily interested in qualitative 
results, the second method with a suitable small time increment h is sufficiently accurate. 

The iteration starts with the initial condition x = x± = . . . = x^-i = 1. Using the 
step size h = 10~ 2 and iterating up to n = 10 5 steps we monitor the distribution of x n . 
As observed in Ref. [10], the x n are positive and fluctuate strongly over several orders 
of magnitude. Since the fluctuations are caused by the multiplicative character of the 
evolution equation, it is reasonable to investigate the distribution of the logarithm 



The distributions for c = 1 and a = 1/2 are shown in Fig. For 6 = the values of y n 
in an ensemble of independent runs are distributed according to a Gaussian. Because 
of the special choice a = 1/2 in the Ito scheme the distribution is centered at the origin. 
As expected, the distribution becomes broader as time proceeds. For b > 0, however, 
the situation is very different. On the one hand the distribution has a smaller width 
and moves to the right as time proceeds. On the other hand it is no longer Gaussian, 
instead it becomes skew (see Fig. [2|). As we will see below, this skewness is one of the 
hallmarks of dynamical scaling in the KPZ universality class. 



y n = log(x„). 



(7) 
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3. Relation to the KPZ equation 

Let us now demonstrate how the stochastic differential equation (J2J) is related to the 
KPZ equation. We start out with the discretized equation ([ED with parameters a, b 
and h: 



ah)x n _i + bhx n ^ k + Vhx n _iz n -i (8) 



Consider the mapping N — > N ® {0, . . . k — 1} : n — > (i, j) 

i = n div k , (9) 
j = n mod k , 

where 'div' and 'mod' denote integer division and modulus. As explained in the 
Introduction, this mapping slices the time series into equidistant blocks of size k. 
The blocks are indexed by i while the index j encodes the position within each block. 
This transformation maps the one-dimensional time series with non-local (time-delayed) 
updates onto a 1+1-dimensional model with local (nearest-neighbor) updates of the 
form [14] 

Xij = (1 + atyx^^ + bhx^j + VAXij-i^j-i , (10) 

where the indices i and j may be interpreted as 'time' and 'space'. Note that the 
boundary conditions differ from periodic ones in so far as they involve a shift in the 
index i: 

Carrying out the continuum limit and performing a suitable Galilei transformation (see 
Appendix) the bulk equation (PLOl ) can be written in terms of differential operators as 

V s x(s, r) = (1 + £)z(s, r) - —rV 2 r x(s, r) (12) 



1 



b' y ' ' lab 



2 

(/ 2a 2 v r 



(x(s,r)£(s,r)) , 



where s and r are temporal and spatial coordinates respectively. Apart from the higher- 
order derivatives in the last term this equation has the same structure as the Langevin 
equation for 1+1-dimensional systems with multiplicative noise [15, 16] which is known 
to exhibit dynamical scaling. This means that the equation is invariant under the scaling 
transformation 

r^Ar, s -> A z s (13) 

with a certain dynamical exponent z > 1. Simple power counting of the leading 
contributions reveals that z = 2 above the upper critical dimension d c = 2. Moreover, 
the higher-order derivatives in the last term of Eq. ( |T2l turn out to be irrelevant in the 
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scaling limit and may be discarded. It has widely believed that this applies also to the 
1+1-dimensional case. The remaining equation reads 

V s x(s, r) = (1 + ^)x(s, r) - ^ V r x ( s > r ) + r )£( s > r ) • ( 14 ) 

With a Hopf-Cole transformation x(s,r) = exp(y(s, r)) this equation turns into the 
stochastic differential equation 



a 1 \ 1 



1 



V s y(s,r) = 1 +\l-2b)-2^ Kv(s,r) + (V r y(s,r)y]+ T £(s,r)(l5) 



b 



with a non- multiplicative noise. Here the additional term — ^ is due to the ltd scheme 
under the transformation. This equation is known as the Kardar-Parisi-Zhang (KPZ) 
equation for non-linear surface growth [11], where y(s, r) denotes the height of a growing 
interface at position r at time s. Since s 

simeqt/r the mapping to the KPZ equation allows us to make the following predictions: 

• The width of a roughening KPZ interface is known to grow as s 1 / 3 , hence the 
variance of y(t) = logx(t) in the original problem should grow as cr 2 (i) ~ t 2 ^ 3 . 

• In the 1+1-dimensional formulation finite-size effects are expected to occur when 
the spatial correlation length £j_ ~ s 1//z becomes comparable with the width of the 
strip r, where z = 3/2 is the dynamical exponent of the KPZ universality class in 
1+1 dimensions. Therefore, finite-size effects are expected to set in when t ~ r 5 / 2 . 

• The (rescaled) distribution of y(t) in the range r <C t <C r 5 / 2 is no longer given by 
a Gaussian, instead the distribution becomes skew. 

The critical exponents as well as the rescaled distribution are universal quantities, i.e., 
they are fully determined by the underlying field theory of the KPZ class. In a recent 
breakthrough, Prahofer and Spohn were able to compute the rescaled height distribution 
analytically [17-20]. 



4. Numerical validation 



In order to verify these predictions, we iterated the discretized version of the original 
stochastic differential equation numerically. Throughout all simulations we used the 
choice b = c = 1. In order to avoid floating-point overflows the parameter a was adjusted 
at the beginning of each simulation by an iterative procedure such that y(t) = \ogx(t) 
does not grow on average. We produced three data sets for step sizes h = 0.02, 0.01, and 
0.005. On the time scales of interest these data sets coincide within statistical errors, 
meaning that deviations due to the discretization can be neglected. 

The left panel of Fig. [3] shows the variance cr 2 (t) = ((y(t) — (y(t))) 2 ) as a function 
of t for two different delays r = 25 and r = 100. As expected, the variance grows faster 
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Figure 3. Variance of y(t) = log(x(t)) in an ensemble of independent iterations 
for b = c = 1 and two different values of the delay r (see text). The dashed line 
indicates the slope 2/3. 




Figure 4. Left: Skewness of the probability distribution of y(t) as a function of 
time for b = c = 1 and r = 100. The theoretical prediction for KPZ roughening 
is marked as a dashed line. Right: Form of the distribution at time t = 30000 
for the same parameters compared to a high-precision profile of a KPZ growth 
process. To demonstrate the skewness we fitted a Gaussian (dashed line) in the 
center. 



for small r. The right panel shows the same data in a double-logarithmic representation 
as a function of the quotient s ~ t/r. For small t/r the two curves collapse and one 
can clearly see the block structure in form of a staircase. The curve for the large 
delay r = 100 displays a clean power law with the slope 0.665(5), confirming the KPZ 
roughening-exponent 2/3 = 2a/ z = 2/3 in 1+1 dimensions. The curve for the smaller 
delay r = 25 begins to deviate from this power law at t/r ~ 100. These deviations are 
caused by finite-size effects and are expected to set in at a time scale t ~ t 5 / 2 . 

In order to substantiate these findings even further, we monitored the skewness 

q(f) to® Ml ~ M*)» 3 > 

1 ' ' °(tf " ((y(f) - (y(tW) 3/2 [ ' 

as a function of time, using the parameters b = c = 1 and r = 100. A Gaussian 
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distribution would yield the value S = while a KPZ profile is characterized by a non- 
zero value. In more recent simulations [21] this value was estimated by S KPZ = 0.28(4) 
while Prahofer and Spohn find the analytical result [17-20] 

S KPZ sa 0.2935. (17) 

As shown in the left panel of Fig. [4[ the skewness of the distribution quickly saturates at 
a constant which is compatible with the theoretical prediction (dashed line). The right 
panel shows the actual form of the distribution at time t = 30000 in comparison with a 
(rescaled) high-precision profile from a simulation of the so-called single-step model [22] 
which is known to belong to the KPZ universality class. As can be seen, both curves 
coincide perfectly and deviate significantly from a Gaussian (shown as a dashed line). 

To summarize, there is strong numerical evidence that y(t) = log(x(t)) is indeed 
distributed in the same way as the heights of a roughening KPZ interface, supporting 
the analytical derivation in Sect. 3. 

5. Conclusions 

In this paper we have shown how a one-dimensional system with time-delayed dynamics 
governed by Eq. ((2j) can be mapped onto a 1+1 dimensional stochastic model with 
local dynamics in the KPZ universality class. Our analytical results are validated by 
numerical simulations. The example discussed here involves two delays. One of them 
is local and arises from the discretization of the time derivative, while the other one is 
nonlocal. One interesting question is whether the scheme here presented can be extended 
to more general systems where both delays are nonlocal 

x{t) =ax{t-Ti)+bx{t-T 2 ) + cx{t)£{t). (18) 

Some recent results show that for such stochastic dynamical rules it is always possible 
to find an equivalent 1+1 dimensional system with local interactions but time-delayed 
boundary conditions, the only condition being that the discretized counterparts of T\ 
and T2 be coprime (two numbers are coprime if they have no common divisor other than 
1). Hence we conjecture that the KPZ universality class of growing interfaces should 
encompass a much broader class of systems whose dynamics is time-delayed. Work in 
this direction is currently under progress [14]. 
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Appendix A. Continuum limit 

In this appendix we show how the continuum limit is carried out. Starting point is the 
1+1-dimensional discretized equation 



x id = (1 + alijXij-! + bh Xi_ 1}j + . (A.l) 

To perform the continuum limit we introduce continuous coordinates s, r by 

s = h s i , r = h r j (A. 2) 

with step sizes h s and h r to be determined later. The discretized field and the noise are 
approximated by smooth functions 



Xi 



„,j -^x(s,r); Zij — > \fh 8 h T £(s,r) , 
where £(r, s) is a Gaussian white noise with the correlations 
(Z{r,sW,s')) = 5(r-r')5(s-s'). 

With this notation the equation above takes the form 

x(s, r) = (1 + ah) x(s, r — h r ) + bh x(s — h s ,r) 
+ a/ hh r h s x(s, r — h r )£(s, r — h r ) , 

This equation may be written equivalently as 

x(s, r) = (1 + ah)e- hrVr x{s, r) + bh e~ hsVs x{s, r) 
+ ^hh r h s e- hrVr [x{s,r)^{s,r)] . 

In order to compensate for a possible drift we apply a Galilei transformation 



r + vs , 



V. s — VVr 



(A.3) 



(A.4) 



(A.5) 



(A.6) 



With this transformation the equation becomes 

e- hsvVr x(s, r) = (1 + ah)e~ {hr+vhs)Vr x(s, r) + bh e~ hsVs x(s, r) 
+ y/ hhrhs e -V* +vh >^ [x(s, r)£(s, r)} . 
Expanding to first order in s and to second order in r we obtain 
bhh s 'V s x(s,r) = 

h(a + b)x(s,r) — (1 + ah)h r + ahh s v V r x(s,r) 

+ - (1 + ah) (hi + 2h r h s v) + ahh 2 s v 2 V 2 r x(s, r) 

_ L 

+ y/ hhrhs (l - (hr + vh s )V r + hyh r + vh s ) 2 V 2 ^j [x(s, r)£(s, r)} . 



(A.7) 

(A.8) 
(A.9) 

(A.10) 
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Now we adjust the velocity v in such a way that the drift (the second term on the rhs) 
vanishes. This can be done by setting 



h r /l + ah\ 
h.\ ah J ' 



leading to the equation 

bhh s V s x(s, r) = h(a + b)x(s, r) — 



2ah 



+ yj hh r h s 



h r 



ah^ r ^ 2 \ah 



In the limit h r = h — > and taking h s — 1 one finally obtains 
V s x{s,r) = (1 + ^)x(s,r) - — ^i(s,r) 



(A.ll) 

(A.12) 

(a;(sr)f(s,r)) 

(A.13) 



+ 



1 



1 + -V r + 



1 

2^2 



V 



(:r(s,r)£(s,r)) 



which is understood to be iterated in the ltd sense. 
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